class: center, middle, inverse, title-slide # sf : LE package pour le traitement et l’analyse de données vectorielles géolocalisées en R - Partie 1 ## Analyse spatiale des agroécosystèmes avec R ### Antoine Casquin ### Inrae, UMR SAS ### 31 Août 2021 ---
<style type="text/css"> .remark-code { font-size: 16px; } .big-code .remark-code { /*Change made here*/ font-size: 150% !important; } .small-code .remark-code { /*Change made here*/ font-size: 65% !important; } /* code and plot side by side */ .left-code { color: #777; width: 38%; height: 92%; float: left; } .right-plot { width: 60%; float: right; padding-left: 1%; } .left-large-code { color: #777; width: 60%; height: 92%; float: left; font-size: 85% !important; } .right-small-plot { width: 38%; float: right; padding-left: 1%; } </style> # Le package sf : A quoi sert-il ? - Réaliser toutes les opérations communes (et moins) que l'on peut généralement éxecuter dans un logiciel de SIG (ArcGIS, QGIS, etc.) sur des données vectorielles spatialisées - Support de tous type de données vectorielles : polygones, lines, polylignes, points etc. <center>  </center> --- # Le package sf : A quoi sert-il ? - Réaliser toutes les opérations communes (et moins) que l'on peut généralement éxecuter dans un logiciel de SIG (ArcGIS, QGIS, etc.) sur des données vectorielles spatialisées - Union, Intersection, Différence, Proximité, Buffer - Projection / Re-Projection - Sélection / Jointures spatiales et par attribut - Exports de couches --- ## Le package sf : comment fonctionne-il ? ```r library(tidyverse) ``` ``` ## -- Attaching packages --------------------------------------- tidyverse 1.3.1 -- ``` ``` ## v ggplot2 3.3.5 v purrr 0.3.4 ## v tibble 3.1.3 v dplyr 1.0.7 ## v tidyr 1.1.3 v stringr 1.4.0 ## v readr 2.0.1 v forcats 0.5.1 ``` ``` ## -- Conflicts ------------------------------------------ tidyverse_conflicts() -- ## x dplyr::filter() masks stats::filter() ## x dplyr::group_rows() masks kableExtra::group_rows() ## x dplyr::lag() masks stats::lag() ``` ```r library(sf) ``` ``` ## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1 ``` --- ## Le package sf : comment fonctionne-il ? ### Basé sur des librairies rapides, open-source, et éprouvées en C/C++ - ** GEOS: Spatial Model and Functions ** - ** Geometries **: Point, LineString, Polygon, MultiPoint, MultiLineString, MultiPolygon, GeometryCollection - ** Predicates **: Intersects, Touches, Disjoint, Crosses, Within, Contains, Overlaps, Equals, Covers - ** Operation **: Union, Distance, Intersection, Symmetric Difference, Convex Hull, Envelope, Buffer, Simplify, Polygon Assembly, Valid, Area, Length, etc. .footnote[ https://trac.osgeo.org/geos] --- ## Le package sf : comment fonctionne-il ? ### Basé sur des librairies rapides, open-source, et éprouvées en C/C++ - ** GDAL: Driver pour une multitude de format de vecteur et de raster ** "**GDAL is a translator library for raster and vector geospatial data formats** that is released under an X/MIT style Open Source License by the Open Source Geospatial Foundation. As a library, it presents a single raster abstract data model and single vector abstract data model to the calling application for **all supported formats**. It also comes with a variety of useful command line **utilities for data translation and processing.** " .footnote[ https://gdal.org/ https://gdal.org/drivers/raster/index.html https://gdal.org/drivers/vector/index.html ] --- ## Le package sf : comment fonctionne-il ? ### Basé sur des librairies rapides, open-source, et éprouvées en C/C++ - ** PROJ: Gestion des projections ** "**PROJ is a generic coordinate transformation software that transforms geospatial coordinates from one coordinate reference system (CRS) to another.** This includes cartographic projections as well as geodetic transformations. PROJ is released under the X/MIT open source license " .footnote[ https://proj.org/ ] --- ## Dépendances de sf **sf** est au coeur de l'écosystème de packages "R Spatial".Auparavant il fallait utiliser rgeos, rgdal, sp etc. ** Comme il réunit et uniformise plusieurs anciens packages, il s'est imposé comme LE package. **  https://r-spatial.org/r/2020/03/17/wkt.html --- ## Charger des données spatialisées dans R A partir d'un fichier(geopackage) ```r reg_bretagne <- st_read("raw_data/shapefile/reg_bretagne.gpkg") ``` ``` ## Reading layer `reg_bretagne' from data source ## `C:\Users\acasquin\Documents\formation\slide\test\raw_data\shapefile\reg_bretagne.gpkg' ## using driver `GPKG' ## Simple feature collection with 4 features and 0 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: -5.142285 ymin: 47.27775 xmax: -1.01564 ymax: 48.9003 ## Geodetic CRS: WGS 84 ``` --- ## Charger des données spatialisées dans R A partir d'un fichier(shapefile) ```r uts_st_malo <- st_read("raw_data/shapefile/Stmalo_25000.shp") ``` ``` ## Reading layer `Stmalo_25000' from data source ## `C:\Users\acasquin\Documents\formation\slide\test\raw_data\shapefile\Stmalo_25000.shp' ## using driver `ESRI Shapefile' ## Simple feature collection with 1291 features and 6 fields ## Geometry type: POLYGON ## Dimension: XY ## Bounding box: xmin: 329505.3 ymin: 6844185 xmax: 344211.3 ymax: 6856797 ## Projected CRS: RGF93 / Lambert-93 ``` - Localisation de la couche - Driver utilisé - Nombre d'objets spatiaux et champs dans la table d'attribut - Type de géométrie - Nombre de dimension - Emprise spatial (bounding box) - CRS (Projection) --- ## Structure d'un objet sf  - Un simple "dataframe" avec une colonne de "géométrie" et des objects inclus qui contiennent des informations nécessaires à la spatialisation. - Il est possible de travailler sur ces données spatialisées comme on travaille sur un simple "dataframe". Une partie de la "puissance" du package sf est son intégration avec le "tidyverse" .footnote[https://r-spatial.github.io/sf/articles/sf1.html#how-simple-features-in-r-are-organized-1] --- ## Plus d'informations sur la projection ```r st_crs(reg_bretagne) ``` ``` ## Coordinate Reference System: ## User input: WGS 84 ## wkt: ## GEOGCRS["WGS 84", ## DATUM["World Geodetic System 1984", ## ELLIPSOID["WGS 84",6378137,298.257223563, ## LENGTHUNIT["metre",1]]], ## PRIMEM["Greenwich",0, ## ANGLEUNIT["degree",0.0174532925199433]], ## CS[ellipsoidal,2], ## AXIS["geodetic latitude (Lat)",north, ## ORDER[1], ## ANGLEUNIT["degree",0.0174532925199433]], ## AXIS["geodetic longitude (Lon)",east, ## ORDER[2], ## ANGLEUNIT["degree",0.0174532925199433]], ## USAGE[ ## SCOPE["Horizontal component of 3D system."], ## AREA["World."], ## BBOX[-90,-180,90,180]], ## ID["EPSG",4326]] ``` --- ## Plus d'informations sur la projection ```r st_crs(uts_st_malo) ``` ``` ## Coordinate Reference System: ## User input: RGF93 / Lambert-93 ## wkt: ## BOUNDCRS[ ## SOURCECRS[ ## PROJCRS["RGF93 / Lambert-93", ## BASEGEOGCRS["RGF93", ## DATUM["Reseau Geodesique Francais 1993", ## ELLIPSOID["GRS 1980",6378137,298.257222101, ## LENGTHUNIT["metre",1]]], ## PRIMEM["Greenwich",0, ## ANGLEUNIT["degree",0.0174532925199433]], ## ID["EPSG",4171]], ## CONVERSION["Lambert-93", ## METHOD["Lambert Conic Conformal (2SP)", ## ID["EPSG",9802]], ## PARAMETER["Latitude of false origin",46.5, ## ANGLEUNIT["degree",0.0174532925199433], ## ID["EPSG",8821]], ## PARAMETER["Longitude of false origin",3, ## ANGLEUNIT["degree",0.0174532925199433], ## ID["EPSG",8822]], ## PARAMETER["Latitude of 1st standard parallel",49, ## ANGLEUNIT["degree",0.0174532925199433], ## ID["EPSG",8823]], ## PARAMETER["Latitude of 2nd standard parallel",44, ## ANGLEUNIT["degree",0.0174532925199433], ## ID["EPSG",8824]], ## PARAMETER["Easting at false origin",700000, ## LENGTHUNIT["metre",1], ## ID["EPSG",8826]], ## PARAMETER["Northing at false origin",6600000, ## LENGTHUNIT["metre",1], ## ID["EPSG",8827]]], ## CS[Cartesian,2], ## AXIS["easting (X)",east, ## ORDER[1], ## LENGTHUNIT["metre",1]], ## AXIS["northing (Y)",north, ## ORDER[2], ## LENGTHUNIT["metre",1]], ## USAGE[ ## SCOPE["Engineering survey, topographic mapping."], ## AREA["France - onshore and offshore, mainland and Corsica."], ## BBOX[41.15,-9.86,51.56,10.38]], ## ID["EPSG",2154]]], ## TARGETCRS[ ## GEOGCRS["WGS 84", ## DATUM["World Geodetic System 1984", ## ELLIPSOID["WGS 84",6378137,298.257223563, ## LENGTHUNIT["metre",1]]], ## PRIMEM["Greenwich",0, ## ANGLEUNIT["degree",0.0174532925199433]], ## CS[ellipsoidal,2], ## AXIS["latitude",north, ## ORDER[1], ## ANGLEUNIT["degree",0.0174532925199433]], ## AXIS["longitude",east, ## ORDER[2], ## ANGLEUNIT["degree",0.0174532925199433]], ## ID["EPSG",4326]]], ## ABRIDGEDTRANSFORMATION["RGF93 to WGS 84 (1)", ## VERSION["EPSG-Fra"], ## METHOD["Geocentric translations (geog2D domain)", ## ID["EPSG",9603]], ## PARAMETER["X-axis translation",0, ## ID["EPSG",8605]], ## PARAMETER["Y-axis translation",0, ## ID["EPSG",8606]], ## PARAMETER["Z-axis translation",0, ## ID["EPSG",8607]], ## USAGE[ ## SCOPE["Approximation assuming equality between plate-fixed static and earth-fixed dynamic CRSs."], ## AREA["France - onshore and offshore, mainland and Corsica."], ## BBOX[41.15,-9.86,51.56,10.38]], ## ID["EPSG",1671], ## REMARK["Parameter values from RGF93 to ETRS89 (1) (code 1591) assuming that ETRS89 is equivalent to WGS 84 within the accuracy of the transformation."]]] ``` --- ## Pourquoi autant d'informations sur la projection ? - Il existe une **multitude** de manières de **représenter** (projeter) la **sphère sur le plan** - **Aucune ne conserve à la fois les angles, les distances et les aires** - Une **projection** est toujours un **compromis**  --- ## Pourquoi autant d'informations sur la projection ? - Le nom de la projection - L'unité (mètre ou degrès) - Un géoïde de référence - Un datum - Un type de projection (conique, cylindrique, etc.) - un code **EPSG** : identifiant unique pour chaque projection --- ## ESPG courants .center[ <table class="table" style="font-size: 20px; margin-left: auto; margin-right: auto;"> <thead> <tr> <th style="text-align:left;"> Nom </th> <th style="text-align:left;"> Unité </th> <th style="text-align:right;"> Code.EPSG </th> <th style="text-align:left;"> Remarque </th> </tr> </thead> <tbody> <tr> <td style="text-align:left;"> World Geodetic System 1984 </td> <td style="text-align:left;"> degré </td> <td style="text-align:right;"> 4326 </td> <td style="text-align:left;"> Coordonnées GPS </td> </tr> <tr> <td style="text-align:left;"> RGF93-Lambert-93 </td> <td style="text-align:left;"> mètre </td> <td style="text-align:right;"> 2154 </td> <td style="text-align:left;"> Projection standard France Métropolitaine </td> </tr> <tr> <td style="text-align:left;"> ETRS89-extended / LAEA Europe </td> <td style="text-align:left;"> mètre </td> <td style="text-align:right;"> 3035 </td> <td style="text-align:left;"> Projection standard UE </td> </tr> <tr> <td style="text-align:left;"> WGS 84 / UTM zone 30N </td> <td style="text-align:left;"> mètre </td> <td style="text-align:right;"> 32630 </td> <td style="text-align:left;"> Universal Transverse Mercator pour longitude entre -6 et 0 </td> </tr> </tbody> </table> ] --- ## Visualiser rapidement un objet "sf" ### En utilisant la méthode "plot" d'un objet sf ```r plot(uts_st_malo) ``` <!-- --> --- ## Visualiser rapidement un objet "sf" ### ggplot ```r ggplot() + * geom_sf(data = uts_st_malo) ``` <!-- --> --- ## Visualiser rapidement un objet "sf" ### Librairie tmap : carte thématique interactive ```r library(tmap) *tmap_mode("view") tm_shape(uts_st_malo) + tm_polygons("PROF") ```
--- class: center, middle # Opérations sur les données vectorielles ## Sélections et opérations sur les attributs ## Opération spatiales : "relations spatiales" --- ## Chargement des données ### BD Carto (IGN), Communes de Bretagne ```r path_bdcarto <-"raw_data/shapefile/BDCARTO_3-2_TOUSTHEMES_SHP_LAMB93_R53_2017-04-18/BDCARTO/1_DONNEES_LIVRAISON_2017-05-00261/BDC_3-2_SHP_LAMB93_R53-ED171/" communes_bretagne <- st_read(paste0(path_bdcarto, "ADMINISTRATIF/COMMUNE.SHP")) ``` ``` ## Reading layer `COMMUNE' from data source ## `C:\Users\acasquin\Documents\formation\slide\test\raw_data\shapefile\BDCARTO_3-2_TOUSTHEMES_SHP_LAMB93_R53_2017-04-18\BDCARTO\1_DONNEES_LIVRAISON_2017-05-00261\BDC_3-2_SHP_LAMB93_R53-ED171\ADMINISTRATIF\COMMUNE.SHP' ## using driver `ESRI Shapefile' ## Simple feature collection with 1233 features and 13 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: 99038 ymin: 6704170 xmax: 400655 ymax: 6885985 ## Projected CRS: RGF93 / Lambert-93 ``` ```r summary(communes_bretagne) ``` ``` ## ID NOM_COMM INSEE_COMM STATUT ## Length:1233 Length:1233 Length:1233 Length:1233 ## Class :character Class :character Class :character Class :character ## Mode :character Mode :character Mode :character Mode :character ## ## ## ## X_COMMUNE Y_COMMUNE SUPERFICIE POPULATION ## Min. :102863 Min. :6706793 Min. : 43 Min. : 67 ## 1st Qu.:219862 1st Qu.:6779488 1st Qu.: 1013 1st Qu.: 690 ## Median :269680 Median :6810069 Median : 1774 Median : 1298 ## Mean :270949 Mean :6807451 Mean : 2207 Mean : 2657 ## 3rd Qu.:330460 3rd Qu.:6837757 3rd Qu.: 2923 3rd Qu.: 2537 ## Max. :398094 Max. :6878456 Max. :16323 Max. :213454 ## INSEE_ARR NOM_DEPT INSEE_DEPT NOM_REGION ## Length:1233 Length:1233 Length:1233 Length:1233 ## Class :character Class :character Class :character Class :character ## Mode :character Mode :character Mode :character Mode :character ## ## ## ## INSEE_REG geometry ## Length:1233 MULTIPOLYGON :1233 ## Class :character epsg:2154 : 0 ## Mode :character +proj=lcc ...: 0 ## ## ## ``` --- ## Sélection des 10 premiers "sf" (lignes) ```r communes_10 <- communes_bretagne[1:10,] plot(communes_10) ``` ``` ## Warning: plotting the first 9 out of 13 attributes; use max.plot = 13 to plot ## all ``` <!-- --> --- ## Sélection d'un attribut (colonne) ```r sf_pop1 <- communes_10[8] # est un objet sf (8ème colone = POP) sf_pop2 <- communes_10["POPULATION"] # est le même objet sf que précédément *sf_pop3 <- communes_10$POPULATION # /!\ est un vecteur (contient juste les données numérique de POP) ``` --- ## Sélection d'un attribut (colonne) .pull-left[ ```r plot(sf_pop1) plot(sf_pop2) ``` <!-- --> ] .pull-right[ ```r plot(sf_pop3) ``` <!-- --> ] --- ## Sélection d'attributs multiples ```r # Objet sf (8ème colonne = POP, 7ème=SUPERFICIE) sf_pop_super1 <- communes_10[c(8,7)] # Même ojet sf que précédément sf_pop_super2 <- communes_10[c("POPULATION", "SUPERFICIE")] # /!\ Vecteur (contient juste les données numérique de POP et SUPERFICIE) sf_pop_super3 <- cbind.data.frame(communes_10$POPULATION, communes_10$SUPERFICIE) ``` --- ## Sélection PAR attribut (tidyverse style) ```r communes_ille_10khab <- communes_bretagne %>% # Sélectionner les communes d'Ille et Vilaine de plus de 10k habitants * dplyr::filter(INSEE_DEPT == 35, POPULATION > 10000) %>% #Ordonner la table par ordre décroissant de population dplyr::arrange(-POPULATION) %>% # Ne garder que les variables qui nous intéressent (geometry non selectionnée) dplyr::select(NOM_COMM, POPULATION, SUPERFICIE) # La géométrie est "collante" elle reste attachée au df, # sauf si explicitement retirée avec st_drop() communes_ille_10khab ``` ``` ## Simple feature collection with 10 features and 3 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: 326573 ymin: 6775013 xmax: 391836 ymax: 6854491 ## Projected CRS: RGF93 / Lambert-93 ## NOM_COMM POPULATION SUPERFICIE ## 1 RENNES 213454 5039 ## 2 SAINT-MALO 45980 3658 ## 3 FOUGERES 20189 1046 ## 4 VITRE 17571 3703 ## 5 BRUZ 17372 2995 ## 6 CESSON-SEVIGNE 17233 3214 ## 7 SAINT-JACQUES-DE-LA-LANDE 12303 1183 ## 8 PACE 11288 3494 ## 9 BETTON 10692 2673 ## 10 CHANTEPIE 10576 1198 ## geometry ## 1 MULTIPOLYGON (((350246 6785... ## 2 MULTIPOLYGON (((331819 6845... ## 3 MULTIPOLYGON (((387509 6811... ## 4 MULTIPOLYGON (((384056 6789... ## 5 MULTIPOLYGON (((345585 6775... ## 6 MULTIPOLYGON (((356273 6793... ## 7 MULTIPOLYGON (((345534 6784... ## 8 MULTIPOLYGON (((341544 6792... ## 9 MULTIPOLYGON (((351092 6797... ## 10 MULTIPOLYGON (((354573 6784... ``` --- ## Rajout d'une colonne, et calcul sur les attributs ### C'est tout simple ! .left-large-code[ ```r # Style "R base" *communes_ille_10khab$SUPERFICIE_KM2 <- communes_ille_10khab$SUPERFICIE / 100 # Style "tidyverse" communes_ille_10khab <- communes_ille_10khab %>% * mutate(DENSITE = POPULATION/SUPERFICIE_KM2) ``` ] .right-small-plot[ ```r plot(communes_ille_10khab["DENSITE"]) ``` <!-- --> ] --- ## Opération spatiales : "geometric confirmation" ### Une famille de méthodes : st_contains, st_crosses, st_disjoints, st_intersects etc. <center>  </center> - Renvoi une liste d'indices pour lesquels la condition est vérifiée (Ne modifie pas les objets) - Voir "cheatsheet" et doc pour l'ensemble de méthodes .footnote[https://en.wikipedia.org/wiki/DE-9IM] --- ## Exemple : Quelles sont les communes voisines de Rennes ? ```r rennes_sf <- communes_bretagne %>% dplyr::filter(NOM_COMM == "RENNES") *voisines_rennes_id <- st_touches(rennes_sf, communes_bretagne) voisines_rennes_id ``` ``` ## Sparse geometry binary predicate list of length 1, where the predicate ## was `touches' ## 1: 812, 899, 911, 914, 920, 921, 924, 926, 937, 946 ``` ```r voisines_rennes_sf <- communes_bretagne[voisines_rennes_id[[1]],] #list voisines_rennes_sf ``` ``` ## Simple feature collection with 10 features and 13 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: 340445 ymin: 6779451 xmax: 361584 ymax: 6800140 ## Projected CRS: RGF93 / Lambert-93 ## ID NOM_COMM INSEE_COMM ## 812 BDCSURCO0000000009594601 PACE 35210 ## 899 BDCSURCO0000000009594939 LE RHEU 35240 ## 911 BDCSURCO0000000009594602 VEZIN-LE-COQUET 35353 ## 914 BDCSURCO0000000009594599 SAINT-GREGOIRE 35278 ## 920 BDCSURCO0000000009594935 CESSON-SEVIGNE 35051 ## 921 BDCSURCO0000000009594236 BETTON 35024 ## 924 BDCSURCO0000000009594936 CHANTEPIE 35055 ## 926 BDCSURCO0000000009594600 MONTGERMONT 35189 ## 937 BDCSURCO0000000009595238 NOYAL-CHATILLON-SUR-SEICHE 35206 ## 946 BDCSURCO0000000009594938 SAINT-JACQUES-DE-LA-LANDE 35281 ## STATUT X_COMMUNE Y_COMMUNE SUPERFICIE POPULATION INSEE_ARR ## 812 Commune simple 345190 6794589 3494 11288 3 ## 899 Commune simple 344426 6788181 1889 8159 3 ## 911 Commune simple 346752 6790110 786 5153 3 ## 914 Commune simple 352140 6794199 1730 9195 3 ## 920 Commune simple 358059 6790333 3214 17233 3 ## 921 Commune simple 354811 6796834 2673 10692 3 ## 924 Commune simple 357489 6785835 1198 10576 3 ## 926 Commune simple 349277 6794313 467 3267 3 ## 937 Commune simple 352921 6782505 2651 6910 3 ## 946 Commune simple 348212 6785736 1183 12303 3 ## NOM_DEPT INSEE_DEPT NOM_REGION INSEE_REG ## 812 ILLE-ET-VILAINE 35 BRETAGNE 53 ## 899 ILLE-ET-VILAINE 35 BRETAGNE 53 ## 911 ILLE-ET-VILAINE 35 BRETAGNE 53 ## 914 ILLE-ET-VILAINE 35 BRETAGNE 53 ## 920 ILLE-ET-VILAINE 35 BRETAGNE 53 ## 921 ILLE-ET-VILAINE 35 BRETAGNE 53 ## 924 ILLE-ET-VILAINE 35 BRETAGNE 53 ## 926 ILLE-ET-VILAINE 35 BRETAGNE 53 ## 937 ILLE-ET-VILAINE 35 BRETAGNE 53 ## 946 ILLE-ET-VILAINE 35 BRETAGNE 53 ## geometry ## 812 MULTIPOLYGON (((341544 6792... ## 899 MULTIPOLYGON (((345728 6785... ## 911 MULTIPOLYGON (((344430 6790... ## 914 MULTIPOLYGON (((349880 6794... ## 920 MULTIPOLYGON (((356273 6793... ## 921 MULTIPOLYGON (((351092 6797... ## 924 MULTIPOLYGON (((354573 6784... ## 926 MULTIPOLYGON (((348171 6795... ## 937 MULTIPOLYGON (((350571 6780... ## 946 MULTIPOLYGON (((345534 6784... ``` --- ---count: false ## Visualisation (ggplot) .panel1-ggplot-voisin-rennes-auto[ ```r *ggplot() ``` ] .panel2-ggplot-voisin-rennes-auto[ <!-- --> ] --- count: false ## Visualisation (ggplot) .panel1-ggplot-voisin-rennes-auto[ ```r ggplot() + * geom_sf(data = rennes_sf, * fill = "green") ``` ] .panel2-ggplot-voisin-rennes-auto[ <!-- --> ] --- count: false ## Visualisation (ggplot) .panel1-ggplot-voisin-rennes-auto[ ```r ggplot() + geom_sf(data = rennes_sf, fill = "green") + * geom_sf(data = voisines_rennes_sf, * fill = "blue") ``` ] .panel2-ggplot-voisin-rennes-auto[ <!-- --> ] <style> .panel1-ggplot-voisin-rennes-auto { color: black; width: 38.6060606060606%; hight: 32%; float: left; padding-left: 1%; font-size: 80% } .panel2-ggplot-voisin-rennes-auto { color: black; width: 59.3939393939394%; hight: 32%; float: left; padding-left: 1%; font-size: 80% } .panel3-ggplot-voisin-rennes-auto { color: black; width: NA%; hight: 33%; float: left; padding-left: 1%; font-size: 80% } </style> .left-large-code[ ```r ggplot() + geom_sf(data = rennes_sf, fill = "green") + geom_sf(data = voisines_rennes_sf, fill = "blue") ``` ] .right-small-plot[  ] --- ## Exemple : Par quelles communes passe la Vilaine ? ```r hydro_bretagne <- st_read(paste0(path_bdcarto, "HYDROGRAPHIE/TRONCON_HYDROGRAPHIQUE.SHP")) ``` ``` ## Reading layer `TRONCON_HYDROGRAPHIQUE' from data source ## `C:\Users\acasquin\Documents\formation\slide\test\raw_data\shapefile\BDCARTO_3-2_TOUSTHEMES_SHP_LAMB93_R53_2017-04-18\BDCARTO\1_DONNEES_LIVRAISON_2017-05-00261\BDC_3-2_SHP_LAMB93_R53-ED171\HYDROGRAPHIE\TRONCON_HYDROGRAPHIQUE.SHP' ## using driver `ESRI Shapefile' ## Simple feature collection with 26171 features and 9 fields ## Geometry type: LINESTRING ## Dimension: XY ## Bounding box: xmin: 102391 ymin: 6705160 xmax: 400947 ymax: 6881302 ## Projected CRS: RGF93 / Lambert-93 ``` ```r # Sélection par attribut (TOPONYME) cours_vilaine <- hydro_bretagne %>% filter(TOPONYME == "fleuve la vilaine") %>% st_union() # Quelles communes (ID) sont "intersectées" par la Vilaine ? communes_vilaine_id <- st_intersects(cours_vilaine, communes_bretagne) # Sélection par ID dans le sf "communes_bretagne' communes_vilaine_sf <- communes_bretagne[communes_vilaine_id[[1]], ] ``` --- ---count: false ## Visualisation (ggplot) .panel1-ggplot-vilaine-auto[ ```r *ggplot() ``` ] .panel2-ggplot-vilaine-auto[ <!-- --> ] --- count: false ## Visualisation (ggplot) .panel1-ggplot-vilaine-auto[ ```r ggplot() + * geom_sf(data = communes_bretagne) ``` ] .panel2-ggplot-vilaine-auto[ <!-- --> ] --- count: false ## Visualisation (ggplot) .panel1-ggplot-vilaine-auto[ ```r ggplot() + geom_sf(data = communes_bretagne) + * geom_sf(data = communes_vilaine_sf, * fill = "cadetblue1") ``` ] .panel2-ggplot-vilaine-auto[ <!-- --> ] --- count: false ## Visualisation (ggplot) .panel1-ggplot-vilaine-auto[ ```r ggplot() + geom_sf(data = communes_bretagne) + geom_sf(data = communes_vilaine_sf, fill = "cadetblue1") + * geom_sf(data = cours_vilaine, * color = "darkblue") ``` ] .panel2-ggplot-vilaine-auto[ <!-- --> ] <style> .panel1-ggplot-vilaine-auto { color: black; width: 38.6060606060606%; hight: 32%; float: left; padding-left: 1%; font-size: 80% } .panel2-ggplot-vilaine-auto { color: black; width: 59.3939393939394%; hight: 32%; float: left; padding-left: 1%; font-size: 80% } .panel3-ggplot-vilaine-auto { color: black; width: NA%; hight: 33%; float: left; padding-left: 1%; font-size: 80% } </style> .left-large-code[ ```r ggplot() + geom_sf(data = communes_bretagne) + geom_sf(data = communes_vilaine_sf, fill = "cadetblue1") + geom_sf(data = cours_vilaine, color = "darkblue") ``` ] .right-small-plot[  ] --- ## Exporter une couche de données ```r st_write(cours_vilaine, "4_output/cours_vilaine.gml") st_write(communes_vilaine_sf, "4_output/communes_vilaine.shp") ``` --- class: center, middle # Exercices ## sf, partie 1 --- ## Voir dans formation/exercice/4_1_sf_Enonce ### Sélection par attribut et relations spatiales simples 1. Créer un objet sf **commmunes_finistere**, qui contient les communes du Finistère (sélection par attribut) 1. Calculer la densité de population (hab/km2) 1. Visualiser (ggplot) la densité de population, exporter le graphe en 300 dpi 1. Exporter la couche sous format `.gpkg` 1. Quelles sont les **3 communes les plus peuplées** du Finistère ? 1. Les afficher en **rouge** sur la carte des communes du Finistère 1. Trouver les communes qui **NE touchent PAS** ces 3 communes 1. Les afficher en **vert** sur la carte 1. Afficher les **communes restantes** en **jaune** 1. Exporter la carte : 160 mm largeur, 120 mm hauteur, pdf --- ## Voir dans formation/exercice/4_2_sf_Enonce ### Imbrications multiples 1. La couche **BV.json** contient les principaux bassins versants bretons. Visualiser ces bassins versants (`fill = red`) en ajoutant `alpha = 0.2` (paramètre de transparence) 1. L'ojectif final est d'obtenir et de représenter deux objets sf qui contiennent seulement des bassins versants indépendants. Le premier objet contiendra la liste des bassins indépendants les plus petits possibles, le deuxième celui des bassin les plus grands possibles. 1. Obtenir la liste des BVs qui ne sont imbriqués dans aucun BV (astuce : regarder la taille de l'élément) 1. Obtenir la liste des BV qui n'incluent aucun autre BV 1. Trouver un exemple de chaque cas suivant : - BV indépendant (n' inclus, ni n'est inclus dans aucun BV) - Tête de bassin versant (n'inclus aucun BV, mais est inclus dans 1 ou plusieurs BV) - BV Intermédiaire (inclus et est inclus dans un ou plusieurs BV) - Grand BV (inclus un ou plusieurs BV, mais n'est inclus dans aucun) 1. Créer, puis visualiser l'objet **petit_bv_indep** qui contient les plus petits BVs indépendants 1. Créer, puis visualiser l'objet **grands_bv_indep** qui contient plus petits BVs indépendants